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ABSTRACT 

Evolution and disruption of galaxies orbiting in the gravitational field of a 
larger cluster galaxy are driven by three coupled mechanisms: 1) the heating 
due to its time dependent motion in the primary; 2) mass loss due to the 
tidal strain field; and 3) orbital decay. Previous work demonstrated that tidal 
heating is effective well inside the impulse approximation limit. Not only does 
the overall energy increase over previous predictions, but the work is done deep 
inside the secondary galaxy, e.g. at or inside the half mass radius in most 
cases. Here, these ideas applied to cannibalization of elliptical galaxies with 
fundamental-plane parameters. 

In summary, satellites which can fall to the center of a cluster giant by 
dynamical friction are evaporated by internal heating by the time they reach the 
center. This suggests that true merger-produced multiple nuclei giants should 
be rare. Specifically, secondaries with mass ratios as small as 1% on any initial 
orbit evaporate and those on eccentric orbits with mass ratios as small as 0.1% 
evolve significantly and nearly evaporate in a galaxian age. Captured satellites 
with mass ratios smaller than roughly 1% have insufficient time to decay to 
the center. After many accretion events, the model predicts that the merged 
system has a profile similar to that of the original primary with a weak increase 
in concentration. 



Subject headings: stellar dynamics — galaxies: kinematics and dynamics - 
galaxies: evolution — galaxies: clusters — galaxies: nuclei — galaxies: elliptical 
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1. Introduction 

The current picture galaxy evolution in clusters leads naturally to galactic cannibalism, 
especially deep in the potential well where the giants reside. Although multiple nuclei 
candidates have been identified (e.g. Tonry 1985, Lauer 1988), recent searches turned up 
many fewer inner core objects than expected (Tremaine 1995). 

A closer look at the evolutionary picture is motivated by the recent demonstration 
that heating of galaxies or star clusters due the time-dependent tidal field can drive 
their evolution at a rate beyond impulse approximation estimates (Weinberg 1994abc, 
hereafter Wl-3 , Murali & Weinberg 1996ab, hereafter MW1-2). This theory invalidates 
the following often-used argument. The higher density of satellite galaxies implies shorter 
internal orbital times than the orbit of the satellite itself. Therefore, the stellar orbits in 
the satellite will adiabatically invariant to the tidal force, and since the dynamical friction 
time scale is much less than a Hubble time for galaxy masses above 1O 9 M , the satellite 
should sink to the center without suffering tidal disruption and remain a distinct compact 
entity. This paper presents estimates of the evolutionary path and evaporation lifetime for 
the cannibalized fundamental-plane elliptical galaxies. We will find that ellipticals with 
sufficient mass to decay are heated and evaporated before a multiple nucleus system can 
result, although such systems may exist transiently. The results also illustrate the interplay 
between tidal heating, tidal stripping and orbital decay. The likely evaporation of accreted 
galaxies may help reconcile the observation of a bimodal velocity distribution of multiple 
nuclei (Tonry 1985) as 1) a dynamical friction mediated selection effect (Merritt 1984) or 
recently captured secondary and 2) a transient population of evaporating secondaries. 

We begin with a description of the astronomical scenario in §|2|. All members of the 
fundamental plane are represented by a spherical model with fixed concentration; this is 
consistent with the observed fundamental plane relations given the Faber- Jackson (1976) 
relation although the best estimates suggest a weak dependence on concentration. We 
assume that the secondary is captured from the cluster by dynamical drag and consider 
evolution after the secondary is bound to the primary. The important dynamical ingredients 
and their implementation is briefly discussed in §^|; the technical details can be found in 
the Appendix and elsewhere (W2, MW2). The results for a number astronomical scenarios 
are described in §[|; these include the survival and evolution as a function of mass, orbital 
decay, and the resulting distribution of stripped stars in the primary. A summary and 
discussion is presented in §||. 



2. Astronomical scenario 
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2.1. Background profile and fundamental plane scaling 

I have chosen a King model for both the primary and secondary. King models with 
logc = 2.35 are representative elliptical profiles (e.g. Mihalas & Binney 1981, Vader & 
Chaboyer 1994) although King models are not good fits in all casesQ. Nonetheless, the 
mass model parameterizes the range of stellar orbital times, and this range determines 
the overall evolution rate from the resonant heating process to be described below. An 
appropriate concentration ensures that a realistic range of orbital time scales are included. 
The conclusions (§|j) are weakly dependent on the inner profile and other fine details of the 
model. 

The radius and mass concentration chosen according to two fundamental plane 
relations. The first is based on the virial theorem and the Faber- Jackson relation, L oc cr 4 
(Faber & Jackson 1976), which results in the following scaling: 

tfocM 1 / 2 . (1) 

and p oc M -1 / 2 . The concentration parameter is invariant under any fundamental 
plane scaling also assuming the Faber- Jackson relation. Therefore all three King model 
parameters, mass, tidal or maximum radius, and concentration, are fixed for each secondary 
of given mass. The second is based on recent observed fundamental plane relations (e.g. 
Pahre et al. 1995, Faber 1995): 

R oc M - 9 . (2) 

To reduce the overall number of parameters in this study, I have chosen to retain the 
Faber- Jackson relation and constant-concentration models, even though recent fundamental 
relations predict that central density and therefore concentration class scales with mass. 
Changes in concentration predominantly change the inner profile and, as noted, only weakly 
effect the overall evolutionary track of the accreted secondary. 

There are three remaining parameters: orbital energy, orbital eccentricity and 
secondary to primary mass ratio. Orbital evolution is determined using local dynamical 
friction which requires the secondary to be inside the primary (see §3.4j ). The initial orbits 
for the secondaries, then, are chosen to have an energy whose circular orbit encloses the 
99% of the primary mass. In other words, we consider evolution just subsequent to capture. 
Eccentricity is parameterized by the ratio of orbital angular momentum to the maximum 
defined by the energy of the orbit, k = J/J max (E) and five values are chosen: 0.1(0.2)0.9. A 
pure circular (radial) orbit has k — 1.0 (« = 0.0). Because a captured elliptical is likely to 



2 The concentration parameter is denned as logc = log w (R ma x / Rcore) ■ 



-4- 



be on a eccentric orbit, we will emphasize the k — 0.1 case. The model profile and location 
of initial orbits in the model is shown in Figure [I]. Finally, each set of 5 orbits is evolved for 
4 different secondary to primary mass ratios: 10~ 4 , 10~ 3 , 10~ 2 , and 10 -1 . 

Dimensionless units are chosen for the King model such that G — M — 1 and total 
gravitational potential energy W = —1/2. For the Wq = 9.5 King model, R core = 0.5 with 
outer radius R ma x — 7.91 in these units. I will take a fiducial central cluster galaxy to have 
M = 10 14 M Q inside of R ma x = 300 kpc and which sets the time scale quoted in years in §[|. 
This fiducial choice is similar to that for M87 (e.g. Binney & Tremaine 1987, Merritt & 
Tremblay 1993). A different choice simply shifts the quoted time scale by the ratio 



T ( A/T \ ^ 2 / R \ 



T V 1O14M 0/ V 300 kpc/ 

For reference, orbital periods for the fiducial scaling whose guiding center radii enclose 
10(20)90% for k = 0.1(0.2)0.9 are described in Figure @. 



(3) 



3. Method overview 

Evolution in the cannibalized ellipticals is caused by the following four interacting 
physical effects: 1) resonant heating and orbit shocking; 2) self-consistent gravity; 3) tidal 
stripping; and 4) dynamical friction. These will be briefly described below and in the 
Appendix. We will see that dependencies in the effect of the four physical processes govern 
the subsequent evolution. 



3.1. Resonant heating 

The orbiting secondary galaxy experiences a differential or tidal force. The combined 
strain and compressive force is time- dependent and can do work on the secondary galaxy. 
If the change in tidal force is rapid compared to internal orbital time scales, a gravitational 
shock, the work can be computed using the impulse approximation. However even if the 
change in strain is slower than internal orbital time scales, significant work may still be 
done: most realistic galaxies will have resonances between the two (or more) internal orbital 
frequencies the external forcing frequency which leads to significant energy and angular 
momentum exchange (Wl-3). 

More picturesquely, the time-dependent force will excite a wake in the secondary. The 
wake will be dominated by a quadrupole or bar-like distortion whose pattern speed is 



0.01 0.1 1 10 
r 



0.01 0.1 1 10 
r 



Fig. 1. — Density, mass and potential for Wo = 9.5 King model in dimensionless units. 
The model is nearly isothermal for 0.05 <^ r <^ 5. The diagram at the lower right shows the 
pericenter and apocenter radii for each orbit (ends of segments) and guiding center (circular 
orbit) radii (solid dots). 
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Fig. 2.— Periods of orbit scaled to a central cluster galaxy with M = 10 14 M and 
Rmax = 300 kpc. Dots represent the orbital periods with the guiding center radius enclosing 
10(20)90% of the primary mass. The values k = 0.1(0.2)0.9 are ordered from bottom to top. 

determined by the external frequency. Similar to torquing by spiral arms, this 'bar' then 
couples to the tidal force, transferring energy and angular momentum to resonant orbits. 
The perturbation-theory-derived heating rates used here are in good agreement with n-body 
simulations (cf. MW2, Johnston et al. 1996) with the advantage of being able to follow a 
weak disturbance without noise. 



3.2. Self-consistent gravity 

By Jeans' theorem (e.g. Binney & Tremaine 1987), an equilibrium of regular orbits 
is described by a phase-space distribution function, f — f(I) where I are the actions (or 
energy and angular momentum for a spherical system). The associated potential and 
density solve the Poisson equation by construction. Although the actions of most orbits are 
invariant to the slowly changing tidal strain, the resonant heating described above changes 
the actions of some small subset of orbits resulting in a slightly out of equilibrium system. 

At regular intervals, the Poisson equation is iteratively solved to maintain equilibrium. 
Because the external force is assumed to do negligible work on an internal orbital time 
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scale, all external perturbations may be temporarily turned off which fixes the actions and 
simplifies solution. So the overal evolution consists of two phases: 1) evolution of phase 
space due to external perturbations in a fixed gravitational potential; and 2) dynamical 
readjustment with all perturbations removed. Practically, a new equilibrium is only 
computed when the changing phase-space distribution implies a 1-2% change to the stellar 
orbits. 

As the equilibrium profile evolves, new orbits become resonant with the external force. 
In this way, a small set of resonant orbits at any one time can change the global structure 
over a number of dynamical time scales. Finally for the results below, the resulting 
equilibrium phase-space distribution is forced to be spherical and isotropic. This is not 
an in-principle demand — the numerical implementation is general — but a choice driven by 
available CPU time. 



3.3. Tidal stripping 

The outer boundary of a secondary is defined by the points at which a star is more 
strongly attracted by the primary. For a circular orbit, this point is the analogous inner 
Lagrange point in the restricted three-body problem. However, for an eccentric orbit, this is 
not an easily parameterized problem; these points change as the secondary orbits resulting 
in foliated stable and unstable regions (e.g. Keenan 1981). N-body simulations suggest that 
setting the boundary to the inner Lagrange point at perigalacticon is a fair prescription. 

The location of the inner Lagrange point scales with the ratio of mean density of the 
secondary to mean density of the primary enclosed with the secondary's orbit. Therefore, 
as the secondary evolves due to time- dependent heating as described in § |3.1| , stars may find 
themselves on the unbound side of the tidal limit. This loss of material also changes the 
equilibrium. If too much material is evaporated, global equilibrium may be lost and the 
smaller galaxy "disrupts" . 

3.4. Dynamical friction 

Finally, the orbit itself is evolving by dynamical friction. For small secondaries, 
Chandrasekhar's dynamical friction formula is an acceptable approximation (Chandrasekhar 
1943, see e.g. Binney and Tremaine 1987). This approximation assumes that the primary is 
infinite and homogeneous with the local value of density and distribution of velocities. The 
drag force is anti-parallel to the motion of the secondary assuming velocity isotropy. For 
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large secondaries, the situation is more complex (e.g. Hernquist & Weinberg 1989, Weinberg 
1989) but the local approximation will be used for simplicity. Because the evolution of 
large secondaries is rapid, it is unlikely that this assumption affects any conclusion. Further 
consequences of the decaying orbit are an increasing resonant heating rate and stronger 
tidal limit, both of which accelerate the evolution. 



4. Evolution of satellite galaxies 

The models and methods of §Q and §[5] are applied to groups of twenty models 
each. Each group of twenty has four mass ratios, M ratio = 10 -1 , 1(T 2 , 1(T 3 , 10" 4 and five 
eccentricities, k = 0.1(0.2)0.9 (cf. Fig. |). The two groups discussed here have guiding 
center orbits which enclose 99% of the primary mass. The first group uses the virial scaling 
and the second uses the observed fundamental plane scaling (cf. §|2|). 

Recall that the physical times quoted below assume a primary mass of 10 14 M inside 
of 300 kpc. Equation ([|) may be used to scale to any desired primary mass and radius. 



4.1. Disruption and Survival 
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Fig. 3. — Contours show remaining mass fraction (key at upper left) as a function of time 
(logarithmic scale) and initial value of k = J/J max - The secondary to primary mass ratio is 
10 _1 . The right hand (left hand) plot shows the virial (observed) fundamental plane scaling. 
The scalloping in the contours is caused by the projection of a finite grid. 



Figures |3|-|6] describe the mass evolution for fundamental plane ellipticals as a 
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Fig. 4. — As in Fig. |3] but for mass ratio 10 
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Fig. 5. — As in Fig. § but for mass ratio 10 3 . Range in time is extended to accommodate 
alternative scalings. 
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Fig. 6. — As in Fig. 
alternative scalings. 



but for mass ratio 10 4 . Range in time is extended to accommodate 
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function of initial eccentricity for the four mass ratios. The contours indicate the mass 
fraction remaining at the time indicated. The vertical axis shows increasing initial orbital 
eccentricity. Galaxies to the right of the 0.05 contour have completely evaporated. Heating 
and stripping is severe for the most eccentric orbits and the highest mass ratios. For the 
ratio lO^ 1 , the orbit decays quickly and has evaporated for all eccentricities by roughly 1 
Gyr (Fig. §). We will see in § |4.2| that disruption occurs near the center of the host galaxy. 

The trends are similar for smaller mass ratios. A 1% secondary (Fig. |]) evaporates 
in 10 Gyr for a nearly circular orbits and in roughly 2 Gyr for an eccentric orbit. A 0.1% 
secondary — a large dwarf galaxy — does not completely evaporate in 10 Gyr even for an 
eccentric orbit, although it is close. Evolution is slower for smaller mass ratios because 1) 
the density of the secondary is larger and therefore couples more weakly to the tidal field; 
and 2) the orbital decay rate, which is proportional to the mass ratio, is slower. 

4.2. Orbital decay 

Figure [7] shows the orbital evolution of initially k = 0.1 orbits for the four mass ratios. 
The decay rate is computed using Chandrasekhar's formula with In A = ln(R max /ri/2) where 
ri/2 is the current half-mass radius of the secondary. Orbital torques are also computed 
in the local approximation with Chandrasekhar's formula. These eccentric orbits becomes 
more circular during its decay (Fig. |||) as previously described by Bontekoe & van Albada 
(1987). Initially, the k = 0.1 orbits with guiding center trajectories enclosing 99% of the 
primary mass have apocenters outside the primary (cf. Fig. |IJ) which may lead to an 
overestimate of the decay time. For the fiducial model (M = 10 14 M Q , R max = 300 kpc), 
only the 10% and 1% mass ratio secondaries can decay into the center in roughly 10 Gyr. 
The decay time for lower mass galaxies is increased by the concurrent mass loss. 

The longer lifetimes for large initial k is correlated with the longer orbital decay times. 
This is shown in Figure |9| which describes the evolution of 1% secondaries for k = 0.3(0.2)0.9 
(the second panel in Fig. |7| is the first member of this sequence with k = 0.1). The decay 
rate is nearly constant for orbits in the inner primary (roughly inside the half-mass radius 
of 60 kpc). The rate increases with decreasing initial eccentricity because the secondary 
spends a larger fraction of its time at higher primary density. This trend decreases the 
spread of decay times with eccentricity, but does not compensate for the slower initial 
evolution of low-eccentricity orbits. 

Similarly, the steep gradient in time across the mass contours in Figures reflects 
the rapid mass evolution which takes place during the final stage of orbital decay. This 
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Fig. 7.— Orbital evolution for initial k = 0.1 orbits for the four mass ratios 
1CT 1 , 1CT 2 , 1CT 3 , 1CT 4 (right to left, top to bottom) with observed fundamental plane scaling 
(cf. Figs. 0— |). Satellites with M ratio = 10" 1 and 1(T 2 have zero mass at the final time 
shown. 
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Fig. 8. — Change in k = J/J ma x for mass ratios indicated for orbits with k = 0.1 initially. 

trend is maintained at the smallest mass ratios although full decay takes longer than a 
galaxian age for the fiducial scaling. 

Combining the results of this and the previous subsection, we reach the conclusion that 
satellites which can fall to the center of a cluster giant by dynamical friction are evaporated 
by internal heating. 

4.3. Distribution of stripped material 

As the secondary is stripped and disrupted, its stars preserve the instantaneous 
orbit and build up the primary as suggested by Richstone (1976). The relative density 
distributions for secondaries on eccentric orbits, k = 0.1, and the four mass ratios are shown 
in Figure [U]. 

The higher-mass secondaries, M rat i = 10 _1 and 10~ 2 , lose mass quickly. Material 
is lost most quickly at pericenter and individual episodes of mass loss during each orbit 
are visible in the outer galaxy. For M ratio = 10 _1 , 90% of the mass is lost within the 
half-mass radius of the primary and 20% within 4 core radii. Overall, the remnant profile 
is steeper than the primary and could be a significant contributor to the inner light after 
a few such events. The two low mass ratio cases, M ratio = 10~ 3 and 10~ 4 , lose mass more 
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Fig. 9. — Orbital evolution for initial k = 0.3, 0.5, 0.7, 0.9 orbits for mass ratios 10 2 , with 
observed fundamental plane scaling. 
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Fig. 10.— Distribution of stripped material from disrupting secondary (solid curve) 
compared with r~ 2 profile (dashed) and background density (dotted) at the four mass 
ratios 10 _1 , 1(T 2 , 10" 3 , 10 -4 (right to left, top to bottom). These use the virial scaling and 
dimensionless units (§ |2.1|) . 
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gradually and the distribution of stars lost in the outer primary is more extended than 
the primary, approximating a r~ 2 distribution. Both have guiding center radii larger the 
primary half-mass radius and have lost roughly 80% and 50% of their total mass at the 
point depicted. 

Overall, these results suggest that mass evaporated from the secondary is distributed 
similarly to and maybe steeper than the profile of the primary. After many accretion events, 
the merged profile will be slightly more concentrated (c.f. Fig. |I0|). 

4.4. Comparision with previous n-body simulation 

Unfortunately, there has been very little reported n-body work on merging for 
self-consistent primary and secondary galaxies with masses differing by orders of magnitude. 
A search of the literature (Quinlan 1996/nociteQuin:96) revealed one similar n-body study 
by Balcells & Quinn (1990) who performed simulations with M Tat i = 0.1,0.2 for rotating 
systems to explore the formation of counterrotating cores. Although the their initial 
conditions and goals are sufficiently different to prevent a precise comparison, the results 
here yield a similar scenario. The orbital decay is rapid with the largest fraction the 
secondary evaporating in or near the core (c.f. Fig. || and the first panel of Figs. [7] and |T0"|.). 
Although Balcells & Quinn conclude that the dense core of the smaller galaxy survives the 
orbital decay and settles in the core, I believe that we describe the same the end result: a 
single merged profile rather than a multiple nucleus system. Clearly, this is subject is ripe 
for additional work. 



5. Summary and discussion 

The major conclusions of this work are as follows: 

• Time-dependent heating will evaporate secondaries with mass ratios as small as 1% on 
any initial orbit well within a galaxian age. The difference with the naive prediction 
that the denser satellite galaxies will invariably survive orbital decay is due to the 
breakdown of the one- dimensional adiabatic invariant in three-dimensional stellar 
systems as described in MW1. 

• Secondaries with mass ratios as small as 0.1% on eccentric orbits are significantly 
evolved and nearly evaporated. Because capture by dynamical drag will preferentially 
produce high-eccentricity companions, this predicts a lower limit: captured secondaries 
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with mass smaller than 0.1% of the primary will not evaporate. 

• Captured satellites with mass ratios smaller than roughly 1% will not decay to the 
center in a galaxian age. This limit may be conservative — even lower mass secondaries 
will decay to the center and evaporate — given evidence gravitational lens results 
(Miralda-Escude 1995) that the optically inferred mass profile of cD galaxies joins 
smoothing on the that of the cluster in general. 

• The profile of the mass lost as the satellite decays is similar to but slightly 
more concentrated than that of the original primary. This implies that that the 
concentration of the cluster giant will gradually increase after many mergers. 

• Altogether, we have the result that satellites which can fall to the center of a cluster 
giant by dynamical friction are evaporated by tidal heating in the process. Disruption 
occurs near the center of the primary; material from both cores combine into a single 
entity. This suggests that true multiple nuclei giants should be rare. This scenario 
does not address the possibility that a massive accretion event will lead to a nuclear 
gas accretion and a burst of star formation (Hernquist & Mihos 1995) and perhaps 
form a second nucleus in situ. 

I thank Sandy Faber, Gerry Quinlan, Chigurupati Murali, Doug Richstone and Scott 
Tremaine for comments and suggestions. This work was supported in part by NASA grant 
NAGW-2224, NAG 5-2873 and the Sloan Foundation. 
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A. Orbit shocking 

Shocking caused by an oscillatory perturbation is a straightforward variant and is 
somewhat easier to compute than a one-shot adiabatic disturbance described in W2. For 
example, if a cluster is dynamically part of the thick disk, then the perturbation has the 
form V p = g{t)z 2 where then g{t) is a periodic function of time. The function g(t) may be 
expanded as a Fourier series in its vertical oscillation period, P: 

oo 

9it) = £ g k e^\ (Al) 

k=— oo 

where u = 2n/P and 

g k = ^ T dtg{t)e-^. (A2) 
i Jo 

The Laplace transform of this Fourier series is trivial and is: 

oo 

^ < A3 > 

fe=— oo 

For physical scenarios (e.g. smooth and continuous mass profiles), g k will converge rapidly 
with increasing \k\. The calculation is analogous for orbit shocking with the following 
changes: 

1. The potential perturbation expansion will be more general than the oc z 2 dependence 
and include all 2nd order moments (all Y 2m terms); 

2. The Fourier expansion of g(t) will have two indices corresponding to the radial and 
azimuthal periods of the cluster orbit. 



See MW2 for details. 

We are only interested in the long-term secular change in the distribution function, after 
any transients have decayed. Following W2, we Fourier-Laplace transform the perturbed 
Boltzmann equation. The secular contribution is second-order distribution function and 
the inverse Laplace transform leads to the desired result. The details of the function g(t) 
are not important which allows one to eliminate it altogether (MW2). Alternatively, one 
can choose a convenient form for g(t), such as a square pulse, and perform the transforms 
explicitly. Either way for time scales large compared to the stellar orbital times, the secular 
change due to heating becomes 

^ = -Ekl 2 i-| i • ^ (ku + 1 • n) 
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where V t \ denotes the action-angle transform of the tidal potential. As described in MW2, 
this expression has the form: 

df df 2 d [ AfT ^df\ 

i = ir = aE{ A(E) ik) (A4) 

which may be solved by standard techniques (e.g. Crank-Nicholson or Chang-Cooper 1970 
schemes) . 

MW3 will describe the effects disk-shocking, orbit shocking and thick-disk shocking on 
the galactic population of globular clusters. 



B. Tidal perturbation 



To compute the effect of an orbit in a galaxy, the Galactic potential may be expanded 
in the cluster frame. The force (inertial) is: 



F 



E 



R+r ~i 

<9 2 $ 



R 



j dxidxj 



x 



(Bl) 
(B2) 



R=R(t) 



where R describes the cluster and r the position of a star relative to the cluster. The tidal 
potential then follows directly: 



1 [ (d 2 § 1 d$\ (R-x) 2 1 d$ 



Vt 2 \ {dR 2 RdR.) R 2 



j r 

RdR 



(B3) 



R=R(t) 



Expanding equation ( [BB] ) in spherical harmonics, perturbed quantities may be computed as 
outlined in §[A[ The non-inertial velocity-dependent forces are not easily incorporated into 
a potential and have been ignored here. 



